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VISCOUS-SHOCK-LAYER SOLUTIONS FOR TURBULENT FLOW OF RADIATING 
GAS MIXTURES IN CHEMICAL EQUILIBRIUM 

By 

E.C. Anderson^ 
and 

James N. Moss^ 

SUMMARY 

The viscous-shock-layer equations for hypersonic laminar and 
turbulent flows of radiating or nonradiating gas mixtures in 
chemical equilibrium are presented for two-dimensional and 
axially-symmetric flow fields. Solutions are obtained using an 
implicit finite-difference scheme and results are presented for 
hypersonic flow over spherically-blunted cone configurations at 
freestream conditions representative of entry into the atmosphere 
cf Venus. These data are compared with solutions obtained using 
other methods of analysis. 
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SYMBOLS 



damping factor [eqs. (15) and (16)] 


C. 

1 

mass fraction of species i, p^/p 

• 

-ft 

mass fraction of element i 

4 


N 


f‘% 

p 

frozen specific heat of mixture, V C.C 

1 p , 1 


^ < 

1=1 


c 

Pf 1 

specific heat of species i, C* ./C*' 

P/i P." 


H 

„2 

defined quantity, 



total enthalpy, H + ^ 



N 


h 

enthalpy of mixture, ^ C.h. 



i=l ^ "■ 



enthalpy of undecomposed ablation material 


h. 

1 

enthalpy of species i, h?/U*^ 

X 


.i 

flow index; 0 for plane flow; 1 for axisymmetric flow 


K 

thermal conductivity of mixture, _C* 

■' ' ^ref p,» 


z 

mixing length [eq. (13)] 


M* 

molecular weight 

• 

' M* 

molecular weight of mixture 

• 


number of species 


1 Le 

1 

j 

Lewis number, p*D^.C*/K’* 

13 p' 


f 

r 

2 

k 



N 


Le,T 


turbulent Lewis nunber 


N 


Pr 


Prandtl number, y*C*/K* 


N 


Pr,T 


turbulent Prandtl number, y*C*/K* 

T p T 


N 


Re 


Reynolds number, 


N 


Sc 


Schmidt number, Ng^ = 


n 


coordinate measured normal to body, 


n 


normal coordinate [eq. (14)] 


P 

Q 


pressure-gradient parameter [eq. ( 17 )) 

2 


pressure, p’ 


(u:)^ 


divergence of the net radiant heat flux, 


net radiant heat flux in n-direction, q*/p*U*^ 




component of radiant flux toward the shock 


- > * 


(-) 


component of radiant flux toward the wall 




c ,v; 


convective heat flux to the wall (eq. (11)] 


radius measured from axis of symmetry to point on body 
surface, r*/r* 


n 


r* 

n 


nose radius 


coordinate' measured along body surface, s*/r* 


temperature, '^*/'^ref 


3 


temperature, (U*) ^/C* ^ 

P t 

free-stream velocity 

velocity component tangent to body surface, u*/U* 
friction velocity [eq. (19) J 

velocity component normal to body surface, v*/U* 

scaled mean velocity component [eq. (18)], v 

shock angle defined in figure 1 

angle defined in figure 1 

normal intermittency factor [eq, (22)] 

boundary-layer thickness 

incompressible displacement thickness [eq. (21)] 

number of atoms of the Ath element in species i 

normalized eddy viscosity, 

eddy viscosity, inner law [eq. (12)] 

eddy viscosity, outer law [eq, (20)] 

transformed n-coordinate, n/n^ 

body angle defined in figure 1 


body curvature 



p 
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p 


o 


molecular viscosity, P*/P*(T*^^) 
eddy viscosity 

coordinate measured along body surface, ^ = s 
density of mixture, P*/P* 

Reynolds number parameter. 


p ^U*r * 

^00 00 |*J 


11/2 


o* Stefan-Bolt n constant 

2 3 quantities defined by equations (4h, 4c, 4d) 
Superscripts : 


j 0 for plane flow; 1 for axisymmetric flow 

quantity divided by its corresponding shock value 
* dimensional quantity 

' total differential or fluctuating component 

" shock-oriented velocity component (see fig. 1) 

Subscripts : 

e boundary- layer edge 

i ith species 


I 


ith el'^r.ient 


s 


shock 


w wall 

® free stream 

values for the solid ablation material at the surface 


INTRODUCTION 

Nuinerica' methods for calculating flow fields with ablation 
products injected into a radiating gas mixture in chemical equi- 
librium have been developed by Sutton (ref. 1) and Moss (ret. 2) 
for the analysis of hypersonic flow over blunt entry probe config- 
urations. The numeri<:a_ solution procedure developed by Sutton 
is applicable to laminar, transitional, and turbulent flow and is 
obtained using a coupled inviscid flow-boundary- layer analysis. 

The inviscid flow solution is determined by an explicit time 
dependent finite-difference scheme similar to the method developed 
by Barnwell (ref. 3) , and the boundary-layer equations are solved 
by use of an integral matrix procedure (BLIMP) developed by Bartlett 
and Kendall (ref. 4). 

Moss' analysis is restricted to laminar flow and solutions are 
determined using an implicit finite-difference scheme developed by 
Davis (ref. 5) for solving the viscous-shock-layer equations. The 
principal advantages of this technique are that the solution is 
direct and that the effects of inviscid-viscous interactions are 
included within a single set of governing equations which are 
uniformly valid throughout the shock layer. 

In the analyses presented by Sutton and Moss, the radi^t-ion 
heat transfer is calculated using the method developed by Nicolet 
(refs. 6 and 7). This radiation model assumes a nongray gas and 
accounts for molecular band, atomic line, and continuum transitions. 

The present report presents the development of a viscous-shock- 
layer analysis applicable to laminar and turbulent flow of radiating 
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or nonradiating gas mixtures in chemical equilibrium. This analysis 
is based upon the viscous-shock-layer analysis applicable to tur- 
bulent flow of perfect gases developed by Anderson and Moss (ref. 8) 
and the laminar viscous-shock layer analysis for equilibrium chem- 
istry developed by Moss (ret. 2). 

Results obtained with the present method of analysis are 
compared with methods which include corrections for inviscid-viscous 
interactions. Solutions are presented for a 120-degree (total 
angle) spherically-blunted cone configuration at freestream condi- 
tions representative of entry into the atmosphere of Venus. 
Heating-rate distributions are compared for a cold wall (freestream 
temperature) nonradiating shock layer and a radiating shock layer 
with injected ablation products. 

The av'ailability of comparative data obtained using methods 
corrected for inviscid-viscous interactions is limited. Consequently, 
the data obtained using the present method of analysis are to be 
considered as preliminary and serve priraarily to establish stability 
of the numerical method. It is emphasized that no attempts have 
beer, made to obtain better agreement with either of the analyses 
used in the comparisons. A more extensive data base is necessary 
to establish the validity of the present solution procedure. 

ANALYSIS 

Governing Equations 

Ihe equations of motion for reacting gas mixtures in chemical 
equilibrium are presented by Bird, Stewart, and Lightfoot (ref. 9) . 

The formulation of these equations in body-oriented coordinates 
appropriate for viscous-shock-layer analysis of laminar flow of 
radiating and nonradiating gases is presented by Moss (ref. 2). 

For turbulent flow, the viscous-shock-layer equations are 
derived using methods analogous to those presentevi by Dorrance 
(r f. 10) for the turbulent-boundary-layer equations and are 
c .pressed in non-dimensional form for the coordinate system shown 
in figure 1 as: 
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Continuity 


3 


(r + n cos 0)^pu + ^ (1 + n*c) (r + n co? 6)^pvj = 


= 0 ( 1 ) 


s- momentum: 


/ u 3u 9u UVK \ 

+ 1 lE = oJ 1 


\1 + nK 3s ’ 3n 1 ■*■ n< / 

1 + nK 3s 

1 3n 


( 2 ) 

UUK / 2k . j cos 5 \ r,,, . +. 3u ll 

~ 1 H n< J \1 + nK r + n cob 0/|_^ ^Tn l + nKj| 


n-momentum: 


. ( y IZ + V ^ \ + 3£ = 0 

^ \1 + HK 9s Dn 1 + nK / 9n 


(3) 


Energy : 


„ ( S-_ 9“ + V - V lE + ,£uiv£_ ^ „2 

^ \1 + nK Ts dn J 9n 1 + nK 


— ((t>i + <t>2 'l»3) 


(rr-^ + _j_juo8_0 \ ^ ^ 

\1 + nK r + n cos G / 


(4a) 


-Q 


where 


<^1 = |1 + e 


Pr 


+ Sr \ ^ 

*^Pr,T / 


(4b) 
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(4c) 


♦ 2 “ 


N 


Pr 


e'*'n 


^Pr " ^ (^Pr,T ■ 


3u _ tju < 
^ 3n 1 + n< 


4> 3 ~ 


Pr 


N 


Le 


■ ^ Xpr^T ("I'®’’' ■ 


N 3C. 

E I’l ^ 

i=l 


(4d) 


and 


li = h + 


(4e) 


Elemental continuity; 


3C 3C 


^ ( 1 + nK 3s ^ ^ ”5rT ) (1 + n<) (r + n 


cos 0)] 


1- X 


3 

■ 5 IT 


(5a) 


(1 + n<) (r + n cos 6)^ 

Sr 


/ t ^Pr \ ‘'l 


■ 


\ _i 


1 9n 


/ 



where 


C •= 


N 

Z 


i=l 


M* 


i\ FF” ^i 


(5b) 


State : 


p = pTR*/M*C*- 


(6) 
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Boundary conditions . The boundary conditions at the shook 
are calculated by usii' , the Ranki'.ie-Hugoniot relations. At the 
wall, the no-slip and no-temperature- jump boundary conditions 
are used; consequently, u^ = 0. The wall temperature and mass 
injection rate are either specified or calculated. For the cal- 
culated mass injection conditions, the ablation process is assumed 
to be quasi-steady and the wall temperature is the sublimation 
temperature of the ablator surface. With these assumptions, the 
expression for the coupled mass injection rate is 


m 


-q 

TT 


* 

c .w 



i=l 




(7) 


For ablation injection, the elemental concentrations at the wall 
arc governed by convection and diffusion as given by the equation 



Precursor effects are neglected while the energy reradiated from 
the surface is included in the radiation transport calculations. The 
net radiative flux, q^, can be represented as the difference of 
two components 


q 


r 




(9) 


At the surface 


q 


( + ) 
r,w 


* 


co*T*** 

w 


( 10 ) 


where e is the emissivity of the ablator. 
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% 

« 


1 

r 

1 

\ 


The heat transferred to the wall due to c >nduction and 
diffusion is 




N 




( 11 ) 


Radiative transpor t. The radiative flux, q^, and the diver- 
gence of the radiative flux, 0, are calculated with the radiative 
transport code RAD, as presented in references 6 and 7. The RAD 
computer code has been incorporated in the present viscous-shock- 
layer computer code (HYVIS) and streamlined for computational 
efficiency . 

The RAD code accounts for the effects of nongray self-absorption 
and radia«-ive cooling. Molecular band, continuum, and atomic line 
transit . • are included. A detailed frequency dependence of the 

absorption coefficients is used for integrating over the radiation 
frequency spcctr'Oiii and the tangent slab approximation is used for 
integrating over physical space. 

Thermodyna m ic and transport prop er ties . The equilibrium 
composition is determined by a free energy minimization calculation 
as developed in reference 11. Thermodynamic properties for specific 
heat, enthalpy, and free energy and transport properties for vis- 
cosity and thermal conductivity are required for each species 
considered. Values for the thermodynamic (refs. 12 and 13) and 
transport properties (ref. 14) are obtained by using polynomial 
curve fits. The mixture viscosity is obtained by using the semi- 
empirical formula of Wilke (ref. 15). 

Eddy-Viscosity Approximations 

A two-layer eddy-viscosity model consisting of an inner law 
based upon Prandcl's mixing-length concept and the Clauser-Klebanof f 
expression (based on refs. 16 and 17) fo ' the outer law is used in 
the present investigation. This model, introduced by Cebeci (ref. 18), 

11 
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assumes that the inner law is applicable for ♦■he flow from the 
wall outward to the location where the eddy viscosity given by the 
inner law is equal to that of the outer law. The outer law is then 
assumed applicable for the remainder of che viscous layer. It is 
noted that the eddy viscosity degenerates to approximately zero in 
the inviscid portion of the shock layer. The degeneracy is expressed 
in terms of the normal intermitcency factor given by Klebanoff 
(ref. 17). The expressions used in the present investigation are 
given in the following sections. 

Inner-eddy-viscosi ty aoprox i mat ion . Prandtl's mixing- length 
concept is stated in non-dimensional variables as 




'6J1 

Cn 


( 12 ) 


The mixing length £, is evaluated by using Van Driest 's proposal 
(ref. 19) stated as 


where 




/ 

.A 

8 = k 1 n 

1 - 

exp I- 
\ 

^)j 


r 


1/2 

♦• np 

n = — 
a\i 

t ’,1 
• w 

L “ 

\'>n 

w_ 



(13) 


(la) 


Here, is the Von K^rindn constant, which is assumed to ha/e a 

value of 0.4, and is a damping factor. 

Cebcci h;ef. 18) suggests that for flows with a pressure gradient, 
the damping factor be expressed as 


f + ■" ^ / 2 

A = 26(1 - 11. 8P) 


(15) 


and for flows with both a pressure gradient and mass injection, 

/ + r 1 

= 26 <- — j^exp (11. 8v^) ” (11.8v^)| 


(16) 
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where 


P 


+ 


-0 


2 





(17) 



(18) 


(19) 


Outer-eddy-viscoaity approximation . For the outer region of 
the viscous layer the eddy viscosity is approximated by the Clauser- 
Klebanoff expression 


■H 







where 




dn 


kj = 0.0168 


and 




/n 



5.5 



( 20 ) 


( 21 ) 


( 22 ) 
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Equation (22) is Cebeci's approximation (ref. 18) of the error- 
function definition presented by Klebanoff (ref. 17). 

For equilibri.*iii flow without radiation, the boundary-layer 
thickness 6 is assumed to be the value of n at the point where 



0.995 


(23) 


and is defined by linear interpolation in an array of local total 
enthalpies. This definition is approximately equivalent to the 
usual boundary-layer definition 

^ = 0.995 (24) 

where u is the local value for the undisturbed inviscid flow 
e 

outside the boundary layer. 

The values of the parameters kj and k 2 in equations (13) 
and (20) depend on the flow conditions being considered, as does 
the constant represented by the value 26 in equations (15) and (16) . 
The values given are used for convenience in developing the numerical 
method . 

For radiating gases, the loss of energy from the shock layer 
makes the total enthalpy definition unsatisfactory. For these 
cases, the boundary- layer thickness is assumed to be that portion 
of the shock layer which contributes 95% of the dissipated energy, 
and is defined by the expression 



M(1 + e'^) 


dn 



M (1 + e'*') ‘ 


dn 


0.95 


(25) 


6 is evaluated numerically within the local iteration loop. 
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This definition shows acceptable agreement with the total enthalpy 
definition when applied to nonradiating flows and should be an 
acceptable definition for radiating gases. 

Method of Solution 

Davis (ref. 5) presented a method for solving the viscous-shock- 
layer equations for stagnation and downstream flow. Moss (ref. 2) 
applied this method of solution to reacting multicomponent mixtures. 
The present meuhod of solution is identical to that of references 
2 and 5. Therefore, only an overview of the solution procedure 
is presented here. 

The numerical computation is simplified by normalizing most 
of the variables with their local shock values. The transformed 
independent and dependent variables are 


n/n 

s 

P = P/Pg 

s 

P = P/Pg 

u/u^ 

T = T/T ) 

s 

v/v^ 

H = H/H^ 

U/Ug 

K . K/K^ , 


(26) 


Since the normal coordinate, n, is normalized with respect to 
the local shock stand-off distance, a constant number of finite- 
difference grid points between the body and shock are used. The 
transiormations relating the differential quantities are 


n ' 

3s “ 8^ " n^ ^ 3n 


where 


n* 


s 


dn 

p 

d? 


-L = _L _L 

3n n^ 3n 


(27a) 


(27b) 

(27c) 
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and 


1 32 
3n" 


(27d) 


32 

3n^ 


The three second-order partial differential equations are linearized 
and written in the standard form for a parabolic equation as 


32w 


0*1 


3W 

9n 


+ 02W + 


“3 


+ 



0 


(28) 


where W represents tangential velocity for the ‘-momentum equation, 
enthalpy for the energy equation, and elemental concentrations for 
the elemental continuity equations. For the energy equation, the 
divergence of the radiative flux is incorporated in the 03 term. 
When the derivatives in equation (28) are converted to finite- 
difference form by using Taylor's series expansions, the resulting 
equations are of the following form: 


A W 
n m,n- 


n u 


(29) 


The subscript r. denotes the grid points along a line normal 
to the body surface, and the subscript m denotes the grid stations 
along the body surface. Equation (29) , along with the boundary 
conditions, constitute a system of the tridiagonal form and can 
be solvel . fficiently. 

A variable grid spacing is used in both the tangential and 

normal directions to the surface so that the grid spacing can be 

made small in the region of large gradients. The order of the 

truncation terms neglected are (first order accurate) and 

either An„An„ , or (An^ - An„ ,). 

n n-l n n-i 

The equations are solved at any body station m in the 
order shown in figure 2. The governing equations are uncoupled 
and the values of the dependent variables are computed one at a 
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time. Each of the second-order differential equations is individually 
integrated numerically by using the tridiagonal formalism (eq. (29)]. 
The global continuity equation is used to obtain both shock stand- 
off distance and the v components of velocity. By integrating 
the global continuity equation between the limits of 0 to 1, a 

quadratic oquatjon for n is obtained. For the v component 

s 

of velocity at n» the global continuity equation is integrated 
with respect to n between the limits ot 0 to n . The pressure, 
p, is determined at station m by integrating the normal momentum 
equation with respect to n between the limits of 1 to n> 

The equation of state is used to determine the density. The solu- 
tion is iterated until convergence is achieved. The solution 
advances to the next body station, m f 1, and uses the previous 
converged solution profiles as initial values for starting the 
solution at station m + 1. This procedure is repeated until a 
solution pass is obtained. 

The first solution pass provides a first approximation to the 

flow field solution because the following assumptions are used. 

2 

The thin shock-layer form of tne n-momentum equation, ^ 

' 1 + ntc 3n' 

IS used; the stagnation streamline solution is independent of 

downstream influence; the term dn /df is set to zero at each 

s 

body station; and the shock angle a is assumed to be the same 
as the body angle 6. These approximations are then removed by 
global iteration. Two solution passes are generally sufficient. 

This solution procedure is programmed for the CDC 6600 computer. 

DISCUSSION AND RESULTS 


Numerical solutions obtained with the present method of analysis 
are compared with an integral boundary-layer solution for a non- 
radiating shock layer with a cold wall boundary condition, and with 
a solution which couples the inviscid flow and boundary-layer 
equations for a radiating shock layer with surface ablation. 

For the nonradiating shock layer with a cold wall boundary 
condition, solutions were determined by Edquist* using the integral 

* The boundary-layer data have not been published and are presented 
by permission of C.T. Edquist, Martin Marietta Corp. , Denver 
Division, Denver, CO. 17 


boundary-layer solution procedure (SHIV) discussed in reference 20. 
Freestream conditions correspond to a typical trajectory point for 
entry into the atmosphere of Venus. The freestream velocity, tempera- 
ture, and density are 10 km/s, 200 K, and 0.01 kg/m^, respectively. 

The atmospheric composition expressed in mole fractions is 0.95 CO 2 
and 0.05 N 2 . The body considered is a 120-degree (total angle) 
spherically-blunted cone having a nose radius of 0.368 m and a 
base radius of 0.66 m. The surface temperature is assumed to be 
equal to that of the freestream. 

The inviscid flow field solution used to specify edge conditions 
for the boundary-layer solutions was determined using a single strip 
integral meth.^d which accounts for the upstream influence of the 
sonic corner. This influence cannot be accounted for in the 
present analysis, but as shown in figure 3, this influence is 
significant only in the region 1.6 £ ^ 

the maximum difference in the surface-pressure distribution com- 
puted using the present method and the single strip integrated 
method is less than 4%. 

Heat-t::ansf er rate distrioutions corresponding to boundary- 
layer solutions for both isentropic expansion edge conditions and for 
edge conditions corrected for vorticity effects are compared with 
the present method of solution in figure 4. It is noted that 
the present solution was obtained assuming instantaneous transition 
from laminar to turbulent flow. The heat- transfer rate correlation 
formula used in the boundary-layer analysis includes a transition 
correction. Both the present analysis and the boundary-layer 
analysis corrected for vorticity effects show a significant 
increase in heat transfer when compared with the boundary-layer 
solution for isentropic expansion edge conditions. The present 
analysis and the corrected boundary-layer analysis differ by as 
much as 30% in the region of fully developed turbulent flow. 
Considering the assumption of local similarity used in the integral 
boundary- layer analysis and the different methods of turbulence 
modeling, the differences between the two methods of solution 
are '-.at excessive. 
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The solution presented for a radiating shock layer with injection 
of ablation products corresponds to the Venus entry conditions used 
in the coupled inviscid flow-boundary-layer analysis presented by 
Sutton (ref. 1). The atmospheric composition expres':ed in mole 
fractions is assumed to be 0.97 CO 2 and 0.03 N 2 . Freestream 
velocity, temperature, and density are 8.8 km/s, 180 K, and 0.0058 
kg/m^, respectively. The geometry considered i'' a 120-degree (total 
angle) spherically-blunted cone having a nose radius of 0.325 m and 
a base radius of C.69 m. The ablator material is carbon-phenolic 
having a composition expressed in mass fractions of 0.11 O, 0.004 N, 
0.851 C, and 0.035 H. 

The surface ablation-rate distribution used in the present 
analysi"? corresponds to that determined by Sutton (ref. 1), and 
is shown in figure 5. Comparisons of the surface-pressure distri- 
butions and shock shapes corresponding to the present analysis and 
that of reference 1 are shown in figures 6 and 7, rei-ijectively . 

The maximum differences in the surface-pressure distributions and 
shock-layer thicknesses determined by the two methods of analysis 
are approximately 3%. Since the inviscid solution does not account 
for displacement effects, the differences noted are expected for 
the specified injection rates. Other properties within the 
essentially inviscid portion of the shock layer show similar 
agreement . 

Comparisons of radiative and convective heating-rate distri- 
butions corresponding to the two methods of analysis are shown 
in figure 8. Differences of 5 to 10% are obtained for the radia- 
tive heating-rate distributions, and convective heating-rate 
distributions differ by 10 to 15% in the region of laminar flow. 

In the turbulent flow region, the agreement between the two methods 
of analysis is unsatisfactory. The reason for the opposite trends 
in the turbulent heating-rate distributions corresponding to the 
two methods of analysis has not been determined. Additional calcu- 
lations will be necessary to resolve these differences. 
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CONCLUDING REMARKS 


The results of the present investigation demonstrate that 
numerically stable solutions to the viscous-shock-layer equations 
can be obtained for turbulent flows of radiating and nonradiating 
gas mixtures in chemical equilibrium. Acceptable agreement between 
the present method of analysis and an integral boundary-layer 
analysis is obtained for a nonradiating shock layer without injec- 
tion of ablation products. The agreement between the present 
method of analysis and a solution which couples the ir.viscid 
flow an'’, boundary-layer equations is unsatisfactory for the case 
of a radiating shock layer with ablation products injected into 
the layer. The limited availability of comparative data obtained 
with other methods of analysis is not sufficient to verify the 
present method of analysis. 
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Figure 1. - Coordinate system. 


Figure 2. 
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Figure 7. - Comparison of shock-layer thickness distribution. 



Figure 8. - Comparison of radiative and convective heating-rate distribution. 
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